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Abstract. A stochastic model is here introduced to investigate the molecular 
mechanisms which trigger the perception of pain. The action of analgesic drug 
compounds is discussed in a dynamical context, where the competition with inactive 
species is explicitly accounted for. Finite size effects inevitably perturb the mean- 
field dynamics: Oscillations in the amount of bound receptors spontaneously manifest, 
driven by the noise which is intrinsic to the system under scrutiny. These effects 
are investigated both numerically, via stochastic simulations and analytically, through 
a large-size expansion. The claim that our findings could provide a consistent 
interpretative framework to explain the emergence of cyclic behaviors in response to 
analgesic treatments, is substantiated. 
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1. Introduction 

Pain in animals, including humans, is triggered by the so-called nociceptors, sensory 
neurons that react to potentially damaging stimulus. Neurology textbooks pQ reports 
on the cascade of successive reactions which are activated by the so called noxious 
stimuli: The peripheral terminals of primary sensory neurons launch the signal, which 
is then transmitted to the spinal and supraspinal nuclei and eventually yields to the 
activation of a matrix of cortical areas that are deputed to the conscious experience of 
pain. 

More specifically, the stimulus originating from a bodily harming menace can be 
directly processed through transduction by the receptors located on the nerve terminals. 
Alternatively, an indirect pathway can take over through the activation of transient 
receptor potentials on keratinocytes or the release of intermediate molecules (such as 
the ATP) which, in turn, act on sensory neurons receptors. In the following we shall 
assume the first scenario to hold, and, though certainly important, disregard other 
mechanisms that might be simultaneously in play. In other words, we simplistically 
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imagine that pain receptors act as effective gates, channeling the route to the involved 
cortical circuits. 

Analgesic drugs relieve the pain by interfering with the peripheral and central 
nervous system. Drug molecules bind in fact their target receptors, and consequently 
inhibit the pain perception. To grasp and visualize the essence of the process, one can 
hypothesize that the bound chemical element occludes the path, by impeding the signal 
transduction through the channel envisioned above. 

Analgesic are commonly used in basic research and clinical practice, but their 
interaction with nociceptory and normal sensory processing remains to be fully 
unraveled. Anesthetics are for instance known to modify the electrical recordings 
measured via evoked potentials (EPs) responses [2] , a powerful diagnostic tools employed 
to monitor and characterize a large variety of central nervous system disorders. EPs are 
elicited by a specific stimulus applied to the e.g. pain receptors and consist in recording 
the induced electrical brain activity, as detected by localized electrodes placed on the 
surface of the head. Furthermore, EPs are also useful in documenting objective response 
to pain [31 H] and can thus prove fundamental to elucidate the molecular processes that 
control anesthetic absorption and metabolization. 

Different analgesic agents have been shown to produce intriguingly distinct effects 
at the level of the EPs [5] . Recorded time series of the solicited electric activity display 
in fact remarkably different patterns, which are generically attributed to the chemical 
specificity of the anesthetic compound. Qualitatively, large, regular, oscillations of 
the electric response manifest, latency and amplitude being peculiar traits, supposedly 
related to the molecular characteristic of the administered drug. 

Furthermore, cycles in the perception of pain have been also reported which might 
be hypothetically driven by similar microscopic processes, the interaction between 
the anaesthetic molecules and their targets playing certainly a role of paramount 
importance. Clearly, the individual experience of pain is also influenced by psychological 
and cultural factors, unfortunately difficult to deconvolve when aiming at resolving the 
objective picture. 

The issue of developing a unique interpretative framework to account for the 
presence of such oscillatory regimes has catalyzed vigorous discussions. The puzzle 
of their existence remains however to be fully understood. 

Current mathematical models [6] approach the problem via deterministic 
paradigms, thus neglecting the crucial role which is certainly played by the noise, 
intrinsic to the phenomenon under scrutiny. These aspects become particularly 
important when accounting for the presence of diverse chemical species, which populate 
the stream flow in a spatially diffusive environment. Different chemical entities may 
compete with the drug molecules and occupy the sites located in close vicinity of the 
receptors, thus effectively hindering the binding event. Under specific conditions, such 
competition sustained by the stochastic component of the dynamics might result in large 
temporal oscillations for the amount of bound receptors, a mechanism which could 
explain the emergence of macroscopic cycles for the sensation of pain in response to 
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medicaments. 

In this paper, we shall speculate on the above scenario by putting forward a network 
of chemical reactions and performing a system-size expansion through the celebrated 
van Kampen theory [7]. This enables us to derive a set of linear equations for the 
fluctuations, with coefficients related to the steady-state concentrations predicted from 
the first-order theory (i.e. the deterministic rate equations). Solutions are identified 
for which the deterministic steady-state occurs via damped oscillations: the inclusion of 
second-order fluctuations leads then to the amplification of sustained oscillations. These 
conclusions are briefly discussed with reference to the existing medical literature. 

2. Description of the model 

Within the simplified scenario depicted above, we shall model the chemical interaction 
between a large, though finite, number of drug molecules (anesthetic), hereby termed 
T, and free receptors Rp which represent their binding target. Following a successful 
binding event, a molecule of the species T disappears, leaving an empty case, hereafter 
labelled E. The population of bound receptor is in turn increased by one unit. 
These assumptions formally translate into the compact chemical notation: 

Rf + T —> Rt + E (1) 

where a stands for the associated reaction rate. The inverse reaction corresponding 
to the spontaneous detachment of the bound component may occur with a certain 
probability which motivates the introduction of the dual relation: 

R T + E R F + T 

The anesthetic molecules T surf in a densely packed environment and certainly 
experience collisions with several other microscopic chemical entities, H, which populate 
the streaming flow. Binary interactions between H and T elements, can occur in 
the close vicinity of the receptors {Rf) location, potentially disturbing and eventually 
interfering with the binding event. As a result of an hypothetical collision, the active 
species T can be ejected by the spatial layer immediately adjacent to the receptor, 
leaving behind an empty case E. Beyond this effect, which stems from purely steric 
interactions, one has to account for possible chemical transformations, which might 
occur when individuals from the H and T species encounter: The active chaser T 

| We here assume that the free T molecule is still chemically active and can thus potentially chase for 
unscreened targets. This working hypothesis can be relaxed leading to conclusions qualitatively similar 
to the ones highlighted below. 

§ In principle it would be extremely useful to dispose of experimental estimates for the reaction rates, 
so to define a realistic range of variability for the free parameters in the model. The most reliable data 
concern the so-called (equilibrium) affinity constant for the case of e.g. the Tramadol, an analgesic 
agent which belongs to the class of synthetic opioid. Depending on the target receptor (and on the 
specificity of the chaser's molecule) the affinity constant is reported to vary of a large amount which 
scans two orders of magnitude (from fraction of unity to hundreds) [8l [9] . 
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H + T -U H + H 
H + T H + E 



Figure 1. The main reaction schemes are depicted. The squares stand for the 
inactive species, while the circles represent the drug molecules. The model is then 
complemented with a set of additional reactions (see equations ([2])— (J5j)) , which accounts 
for the possibility that T and H enter (resp. leave) the region deputed for the 
interaction. 

can lose its ability to bind the target [JJ and it is thus mapped into an inactive H 
molecule. To incorporate these effects into the proposed description we postulate the 
following interaction rules, which are loosely inspired to the predator-prey competition 
mechanism: 

H + T H + H 
H + T H + E 

The values of a and 7 characterize the effectiveness of the interaction, which is in turn 
sensitive to the choice of the compound T. The idealized cartoons of figured] are aimed 
at visualizing the above reaction schemes. 

To complete the model we introduce an effective migration, by requiring that the 
T and H molecules can enter (resp. leave) the region deputed to the interaction. The 

|| Note that this can happen both due to a mechanical stress or via chemical combination of the colliding 
species, see for instance [10] where the plasma protein binding is discussed. For a specific application 
relative to the case of the Tramadol refer to [TT] . 
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latter assumption yields to the following set of chemical relations: 

(2) 
(3) 
(4) 
(5) 

The population, namely the ensemble of elements belonging to an homologous species 
X, will be labelled in the following with the symbol nx- Notice that the number of 
receptors N\ = n Rr + n Rp and the total amount of molecules (including the empties) 
N 2 = ut + nn + ue are conserved quantities. This observation enables us to reduce the 
complexity of the problem by setting: 

n Rp = N 1 - n Rr n E = N 2 - n T - n H 

In the following we shall use the vectorial notation n = (ny, tir t , n#) to help keeping 
the mathematical developments compact. 

We are now in a position to define the transition rates T(n |n) from a state n to 
a different one n. In our convention initial states are on the right and final states 
on the left. As an example, consider equation ([T]). The probability to pick up a 
T constituents follows from simple combinatorics and reads ut/N 2 , while there is a 
probability (Nx — n Rr )/Nx of Rp being chosen. This results in a(n T /N 2 )(Nx — n Rr )/Nx 
for this particular transition rate [12]. A complete listing of transition probability 
associated to the preceding set of chemical reactions is here enclosed: 

n T Nx - n Rr 



T(n T — 1, n Rr + 1, n H \n) = a 
T(n T + 1, n RT - 1, n H \n) = /3 



N 2 ^ 
N 2 -n T - n H n Rr 

N 2 Ni 



rr( 1 , 11 \ n Hn T 

T(n T - l,n RT ,n H + l\n) = Ttti^T 

l\ 2 JV; 



2 

rpi -i I \ n H n T n T 

T(n T -l,n RT ,n H \n) = a —— + 6 1 — 

T(n T +l,n RT ,n H \n) =t)x 

l\ 2 

T(n T ,n RT ,n H - l|n) =^2^7- 

l\ 2 

N 2 — up — nu 

T(n T ,n RT ,n H + l\n) = rj 2 

1\ 2 

Having defined the transition rates, the master equation governing the evolution of 
the discrete stochastic model takes the form: 

-^P(ri, t) = T(n\n T + 1, n Rr - l,n H )P(n T + 1, n Rr - 1, n H , t) 
+ T(n\n T - 1, n RT + 1, n H )P{n T - 1, n Rj , + 1, n H , t) 
+ T(n\n T + 1, n Rr , n H - l)P(n T + 1, n Rr , n H - 1, t) 
+ T(n\n T + 1, n Rr , n H )P(n T + 1, n RT ,n H , t) 
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+ T{n\n T - 1, n RT ,n H )P(n T - 1, n RT ,n H , t) 

+ T(n\n T , n RT ,n H + l)P{n T , n Rr ,n H + 1, t) 

+ T(n\n T ,n RT ,n H - l)P(n T ,n RT ,n H - l,t) 

- T(n T - 1, n KT + 1, 7iH-|n) + T(n T + 1, »r t - 1, 7iff|ri) 

+ T(n T - l,n RT ,n H + l|n) + T(n T - l,n RT ,n H \n) 

+ T(n T + l,n RT ,n H \n) + T(n T ,n RT ,n H - l|n) 

+ T(n T , ni? T , n H + 1 |n) -Pfe 



(6) 



where P(n,t) is the probability to find the system in the state n at time t. In the 
next section we shall shortly report about our stochastic simulations, before turning to 
develop the analytical framework. 



3. Numerical simulations 



Based on the chemical equations introduced above, numerical simulations can be 
carried on, which respect the intrinsic stochastic nature of the model. To this end we 
employ the celebrated Gillespie algorithm [13] which exploits the information encoded 
in the reaction scheme to advance the system in time through a sequence of random 
increment^ A randomly selected reaction is forced to occur during each successive 
step. It should be emphasized that time increments and associated reactions are chosen 
so to recover the exact probability distribution of the stochastic time series. For a more 
detail account on the philosophy of the integration recipe, the reader can consult [T3] . 

A typical result is represented in figure [2] where the normalized population of T- 
like molecules and bound receptors Rt are reported, as function of time. Notice that 
large stochastic oscillations are clearly displayed, despite the relatively large number of 
simulated molecules. Even more interestingly, the oscillations persist when increasing 
the populations amount and only when a very large number of constituents is allowed, 
the system eventually sets down to its mean-field solution. 

As already remarked in [T2l HU [15] , this result is odds with the intuitive believe 
that fluctuations can be safely ignored, due to the usual statistical factor 1 / y/N. 
Indeed, the observed fluctuations arise from the amplification of the intrinsic noise, 
associated to the discrete component of the dynamics and can potentially bear and 
extraordinary conceptual relevance in several applications. With reference to the case 
at hand, emerging, regular, oscillatory patterns in the Rt quota, could potentially 

% The standard implementation of the Gillespie algorithm is based on a nested sequence of operations 
which is here briefly recalled. First, one initializes the system at t = 0, by assigning a number of 
molecules to each of the considered species. Then the following iterative scheme runs: (i) Calculate the 
transition rates Ti(n'\n) associated to each prescribed reaction i. The quantity Tq = Yli=i Ti(n'\n) is 
stored; (ii) Extract two random numbers, respectively r\ and r 2 , from a uniform distribution, which are 
used to (a) update the simulation time by a finite amount St = l/Toln(l/ri) and (b) select the index i 
labelling the next reaction which is deputed to occur (i is such that 2fc=i Tk < r^To ^ Sl=i (iii) 
Update the species accordingly and go back to point (i). 
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Figure 2. Drug molecules (panel (a)) and bound receptor (panel (b)) densities as a 
function of time. The solid lines refer to the time series from stochastic simulations, 
while the dashed lines arc calculated from numerical integration of the mean field 
equations dZ])-©. The insets represent a zoom of the initial evolution and allow to 
appreciate the agreement with the mean-field solution at short times. Parameters 
values for both panels arc a = 0.008, (5 = 0.005, 7 = 0.3, St = 0.001, 6 2 = 0.05, 
r?i = 0.001, a = 0.06, Ni = 5300 and N 2 = 20000. 



explain the modulations in the pain perception as reported in the literature, higher 
values of ur t corresponding, in our interpretative scheme, to a less pronounced pain 
level. Alternatively, on a different time scale, the effects here discussed could provide a 
consistent microscopic picture to understand the presence of quasi-periodic fluctuations 
in the evoked electric activity of laboratory animals under anesthetic treatment. 

In the forthcoming sections, we shall gain some analytical insight into the model 
results and characterize the specific traits of the observed oscillatory behaviors through 
the elegant van Kampen [7j expansion. 

4. On the deterministic limit 

The deterministic counterpart of the governing master equation is straightforwardly 
obtained as follows. Focus on T and observe that, by definition: 

(n T ) = ^ n T P(n, t) 

n 

Multiplying equation (jBJ) by n T and summing over n returns on the right-hand 
side d (ut) / dt. Simplifying the left-hand side is somehow more laborious and requires 
some effort. Proceeding in a completely analogous fashion for n# and Ur t yields to the 
following system of coupled differential equation for the ensemble independent variables: 

d /utNx- u Rt \ I ' N 2 -n T -n H n Rr \ I ' n H n T \ 
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\N 2 / ,x \ N 2 
d = a /nr^i-n^X _ ^ / N 2 - n T - n H n Rl 



dr nTI \N 2 Ni I \ N 2 N x 
d /n H n T \ / N 2 - n T - n H \ , / n H 



dt \N 2 N 2 / ' \ N 2 I \N 2 

The mean-field approximation corresponds to ignore the correlations, in the above 
rate equations when performing the limit for N\ and N 2 large. Introducing <pT = 
{ut)/N 2 , (j)R T = (ur t )/Ni, 4>h = (nn)/N 2 , rescaling time as r = t/N 2 and formally 
sending N\, N 2 — > 00, one finally obtains: 

= _ acj) T {l - (j) RT ) + p(p RT (l - (f) T - (j) H ) - (7 + a)(j)H<j>T (7) 

- 5i(j)T + - 4>t - 4>h) 

-t-4>Rt = c [«0t(1 - 4>R T ) - P^Rri 1 -<fo~ 4>h)} (8) 

dr 

—(f) H = "i<p H <pT + 772(1 -<f>T — 4>h) - S 2 4> H (9) 
dr 

where c = N 2 /N\. 

We shall hereon limit our discussion to the case 7] 2 = which will prove analytically 
tractable. The conclusions here demonstrated with reference to the selected case study, 
will obviously apply to the more general setting where fresh H molecules are allowed 
to enter the interaction region. Investigations on the complete model (t] 2 7^ 0) will be 
object of a forthcoming publication. 

Two fixed points, respectively labelled 0* and 0*, are identified: 



—2 



Vl (XT]! Q 



^1 + Vi + ar li 
S 2 a[5 2 (7 + °) + Vil] Viil - 82) - 



7 ' «[5 2 (7 + er) + 7717] + + cr)(7 - 5 2 ) + ^7] ' 5 2 (7 + cr) + 7717 
In figure [21 the solid lines represent the above equilibrium solution: For such choice 
of the parameter, as previously mentioned, the stochastic dynamic displays macroscopic 
oscillation for the monitored quantities, the average reference value being correctly 
predicted by the mean-field theory. What is the reason for these regular fluctuation 
to manifest? Are they resulting from the intrinsic finite size nature of the simulated 
medium? 

To answer these questions it is crucial to determine the stability matrix associated 
with the fixed points, as it will play a central role in investigating the cycling 
phenomenon. One can thus re-write the mean-field equations in the compact form 
dtfik/dr = fk{4>), where the index k = 1,2,3 codes the different species, namely T 
(k = 1), Rt {k = 2) and H (k = 3). The 3x3 Jacobian matrix ,7^ = dfi/d<f)j\FP (here 
FP means "evaluated at the fixed point") controls the linearized dynamics about the 
fixed point and can be cast in the form: 

-a(l — 0Jj T ) — 0tt>R - (7 + a)<S>* H — Si - m a <t>T + /3(l - 4>t ~ 0ff) -P4>* Rt - (l + <*)4>T 
J(±*)=\ c[a(l-</4 r ) + 00£ T ] -c[^ + |3(l-^ -4>h)1 =»r t (10) 

14>"h ~t<t>T ~ S 2 
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One can easily show that J(<f)*) admits two real negative eigenvalues and a third 
one, also real, whose sign depends on the choice of the parameters. A stability analysis 
for the second equilibrium point proves technically difficult. However, via numerical 
inspections, a large region in the parameters' space is identified, which yields to complex 
solutions. In particular, complex eigenvalues of the J(4> 2 ) are found having negative 
real part. This condition ensures an oscillatory approach to equilibrium, a fundamental 
ingredient which is eventually responsible for the large scale modulations observed in the 
finite size regime. Tracing the region in space deputed to the aforementioned behavior 
falls out of the scope of the present paper and shall be postponed to a forthcoming 
contribution together with a detailed characterization of the general case with r) 2 7^ 0. 

Following the above, from now on, we shall denote with <p* a particular value of 0* 
for which damped oscillations do occur in the mean-field scenario. 



5. Characterizing the fluctuations: The van Kampen expansion 

As clearly depicted in figure [2] the innate discreteness of the stochastic medium drives 
into the system important effects which cannot be captured within the continuous mean- 
field scenario. To shed light into the observed phenomena we can bring into the game 
the fluctuations by performing the following explicit replacement: 

n T = N 2 <p T (t) + y/W 2 £, T n Rr = N l( f> RT {t) + v^Gjt n H = N 2 <p H (t) + ^N^h (H) 

where the new continuous variable £ = (£,t,£r t £r h ) replace the discrete quantity n = 
(nT,nR T ,nn) in the definition the probability distribution, namely P(n,t) = II(£, r). 

Before proceeding, it is worth emphasizing that the 1/y/Ni (resp. l/y/N^) term 
holds because of the central-limit theorem: The fluctuations are in fact expected to decay 
in a similar fashion and, in the continuous limit, N\,N 2 — > 00 the system is entirely 
characterized in term of its continuous variables = (4>t, 4>R T i 4>h) as prescribed by 
equations (TTTT) . 

The master equation can be re-written as function of the new variables. The left- 
hand side reads: 

d „, N 1 dU 1 d , dU c" 1 d , dU 1 d , dU 
-P(n,t) = — — 



dr N 2 dr s/W 2 drd^ T v^T dT ^flr v^d-r 

where use has been made of the fact that the time derivative is taken at constant n. 
The right-hand side follows from a straightforward, though lengthy, application of the 
large-iV van Kampen expansion. The main step of the derivation are reviewed in the 
appendix. The interested reader can refer to [3, [16] for a detailed account on the whole 
procedure. 

At leading order of the expansion we recover the deterministic mean-field equations 
flU) - ©, while at next-to-leading order we obtain the linear multivariate Fokker Planck 
equation (1A.5j) that governs the evolution of the fluctuations. The coefficients of 



this equation are completely specified as function of the model's parameters (see the 
appendix): In principle, by solving equation ( 1 A. 51) we are in a position to quantify the 
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deviation from the ideal mean-field dynamics, via the probability distribution IT. At 
present, we aim at understanding the oscillation and, to this end, we invoke a completely 
equivalent formulation of the Fokker Planck equation. The problem can be in fact cast 
as a set of stochastic differential equations of Langevin type, which take the explicit 
form: 

^ = A^) + Vi(r) * = !,.. ,3 (12) 

where for convenience, as a natural extension of the notation introduce above, we have 
now set £i = £t, £2 = £r t , £3 = £,h and A, is specified in the appendix. The term rji is 
a Gaussian noise with zero mean and with correlation given by 

( Vl (T) Vj (r')) = B lJ d(r-r') 

To highlight the existence of a possible oscillatory behavior we perform a Fourier 
analysis of equations ( 1T21 . and obtain: 

j 

where the tilde stands for the Fourier transform. Following [15] we can re-write this as: 

= (13) 

3 

with = —iuiSij — Mij. In addition one gets: 

(fji{u})fjj(u)} = Bij{2ii)5{u + J) 
Solving equation (fT3l) for £j and computing the power spectrum results in: 

j k 

where we have used = §ji(—uj). The power spectrum predicted by equation (|T4j) 

is plotted in figure El for the same set of parameters as employed in the simulations of 
figure [2j A clear peak is detected. Moreover, the theoretical curve interpolates correctly 
the numerical profile. These results confirm that the macroscopic oscillations which 
manifest in the acquired time-series stem from the noise intrinsic to the system under 
investigation. It is our believe that mechanisms similar to the ones here hypothesized, 
are potentially in place in the complex human (animal) body environment and could in 
principle form the basis of a consistent molecular interpretation for the large collection 
of experimental, biomedical observations to which we made reference in the introductory 
section. 

6. Conclusions 

In this paper we propose a discrete dynamical framework which is aimed at shedding 
new light onto the complex molecular processes intervening in response to an external 
harming stimulus so to trigger the pain sensation. We are in particular interested in 
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Figure 3. A plot of the power spectrum P(w) for the Rt time series as a function 
of the frequency to. The noisy line corresponds to the spectrum calculated from 500 
runs of the stochastic simulation of the model. The smooth line shows the analytic 
prediction from equation (|14[) . For the parameters' setting refer to the caption of figure 

m 

elucidating the crucial interplay between the administered drug molecules, which express 
their analgesic function chasing the target receptors, and other chemical elements freely 
diffusing in the stream. The latter can substantially reduce the anaesthetic effect, 
by hindering the available binding sites. Similarly, drug molecules can be turned 
into inactive species following binary encounters. The mechanisms here postulated 
are formally coded via chemical reactions and define a consistent stochastic scheme. 
Numerical simulations display macroscopic oscillations in the concentration amount: 
the number of bound receptors change cyclically in time, a trend which we assume to 
induce an analogous modulation for the experienced perception of pain. These findings 
are analytically illustrated by developing a large-size expansion which enables us to 
predict the existence of a peaked power spectrum. It is important to remark that 
the amplification process here discussed stems from the underlying stochasticity, which 
excites the resonant frequencies of the system. Oscillations arise hence naturally, driven 
by the noise which is intrinsic to the system and without invoking any ad hoc couplings 
among the molecular agents participating to the dynamics. Our findings could suggest 
the existence of a simple, though general, molecular mechanism responsible for the 
emergence of cyclic behaviors in response to analgesic treatments [5]. We shall here 
stress that our main conclusions apply also for the more general setting where r]2 ^ 0. 
In particular the peaked power spectrum is also found in this latter case and the region in 
the parameters' space which corresponds to the emergence of the cycles can be partially 
identified on the basis of explicit analytic formulae. These findings will be reported in 
a forthcoming publication [17] . 
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The first technical point of the van Kampen expansion concerns the introduction of the 
so-called shift operators, Ej; 1 , E^ 1 which obey: 

E^/Gl t) = f(n T ± l,n RT , n H ) 
E Sjfe t) = f(nr, n RT ± 1, n H ) 
E^/fe *) = f( n T, n RT , n H ± 1) . 
The master equation ([6]) is hence cast in the form: 



1 dU 

W 2 ~dr 



dU 
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dU 
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The advantage of using the shift operators relies in that they admit a simple expansion 
in the limit for N\ (resp. N 2 ) large: 







E± 1 = l±N 2 ~ 1/2 

9& 



-1 V 



8^ 



1/2 



T 2 

Q 1 



± • • 



(V 



(A.2) 
(A.3) 

" 2 dfr"^' 2 afe^"" (A ' 4) 
Plugging (lA~2l) -( fA~4jl into (!A~T|l . after some algebraic manipulation, one recovers at the 
leading order the mean-field equations, formally identical to the ones reported above, 
see equations (EJ-tEJ). The next-to-leading order result in a Fokker Planck equation for 
the probability distribution n(£, t): 



E 



±i 

H 



where: 



dU d lv^ n d 2 Tl 



(A.5) 
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with 

/ -a(l-<fi* RT )-P<t>R T -(l + <r)<i>H-ti-m c 1 ^^* -1-/3(1-0*-^)] -p<fi* RT _( 7 + -)^« _ m \ 
M= c i/2 [a(1 „^ ) + / 3^ T ] -c[a0J + /9(1 - - ^)] c 1/2 W>^ T 

\ i<t>* H 7$*, - 5 2 / 

Notice that for c = 1 (see [33]), M reduces to the Jacobian matrix (fTUj) which can be 
directly calculated from the mean-field system. B is instead a symmetric matrix whose 
elements reads: 

B u = a<t>* T {l - <f)* RT ) + (3<t>* RT (l -<P* T - <P* H ) + (7 + + 5^ + Vl (l - - <j>* H ) 

Bl2 = -c^[a<P* T {\ - <P Rt ) + P^ Rt (1 - 4>t - <Ph)\ 
B 13 = -^<f)* T (j)* H 

B 22 = c[a<P* T (l - <P Rt ) + (3^ Rt {1 -<&- 4>h)\ 
B 23 = 

B33 = + 
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